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Abstract. We investigate the nature of the critical behavior of the random- 
anisotropy Heisenberg model (RAM), which describes a magnetic system with random 
uniaxial single-site anisotropy, such as some amorphous alloys of rare earths and 
transition metals. In particular, we consider the strong-anisotropy limit (SRAM), 
in which the Hamiltonian can be rewritten as the one of an Ising spin-glass model 
with correlated bond disorder. We perform Monte Carlo simulations of the SRAM on 
simple cubic lattices of linear size L, up to L = 30, measuring correlation functions of 
the replica-replica overlap, which is the order parameter at a glass transition. The 
corresponding results show critical behavior and finite-size scaling. They provide 
evidence of a finite-temperature continuous transition with critical exponents rj — 
—0.24(4) and v Q = 2.4(6). These results are close to the corresponding estimates 
that have been obtained in the usual Ising spin-glass model with uncorrelated bond 
disorder, suggesting that the two models belong to the same universality class. We 
also determine the leading correction-to-scaling exponent, finding u> = 1.0(4). 



PACS numbers: 75.50.Lk, 05.70.Jk, 75.40.Mg, 77.80.Bh 

1. Introduction 

The critical behavior of magnetic systems in the presence of quenched disorder has 
been the subject of extensive theoretical and experimental study. An important class of 
systems consists in amorphous alloys of rare earths with aspherical electron distributions 
and transition metals, for instance TbFe2 and YFe2- They are modeled pQ by a 
Heisenberg model with random uniaxial single-site anisotropy defined on a simple cubic 
lattice, or, in short, by the random- anisotropy model (RAM) 

H = -J^s x - s y -D^2(u x ■ 4) 2 , (1) 

(xy) x 

where s x is a three-component spin variable, u x is a unit vector describing the 
local (spatially uncorrelated) random anisotropy, and D the anisotropy strength. In 
amorphous alloys the distribution of u x is usually taken to be isotropic, since, in the 
absence of crystalline order, there is no preferred direction. 

Random anisotropy is a relevant perturbation of the pure Heisenberg model, so that 
random-anisotropy systems show a critical behavior that is different from the Heisenberg 
one. For small D however, Heisenberg behavior may still be observed, the Heisenberg 
fixed point controlling the multicritical behavior of the system. If T p = T C (D = 0) is 
the critical temperature of the pure Heisenberg model, in the limit t p = T/T p — 1 — > 
and D — > the free energy has the scaling form J2] 

F=\t p \ 2 - aH f(D 2 \t p \-^), (2) 

where an = —0.1336(15) is the specific-heat exponent in the pure Heisenberg theory 
[HIE] 4>d = 0.412(3) is the crossover exponent associated with the random-anisotropy 
perturbation [2], and f(x) is a universal function.^ 

J As a consequence of Eq. |J2J, for sufficiently small D the critical-temperature shift is given by 
T C (D) - T c (0) w cD 2 ^ + aiD 2 + a 2 D A + ... where 2/0D w 4.9. Note that the nonanalytic term 
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The critical behavior in the presence of random anisotropy has been investigated 
at length, but a satisfactory picture has not been achieved yet. A recent review 
can be found in [Sj. The Imry-Ma argument [01 Ej forbids the appearance of a 
low-temperature phase with nonvanishing magnetization for d < 4. This is also 
supported by field-theoretical renormalization-group (RG) calculations using the replica 
method [5J E3 EH I2|- However, a glassy transition with a low-temperature phase 
characterized by quasi- long-range order (QLRO), i.e., a phase in which correlation 
functions decay algebraically, is still possible [Zj . A Landau-Ginzburg calculation ^T] of 
the equation of state for D — > and a recent 4 — e study |T2J HJ| based on a functional 
RG approach support this scenario in the weak-anisotropy limit. In the large- anisotropy 
limit D — ► oo the model becomes an Ising spin glass with a correlated bond distribution. 
An interesting hypothesis, originally put forward in |14j . is that in this limit the RAM 
transition is in the same universality class as that of the Ising spin-glass model (ISGM) 

[13 CHI El. 

Numerical simulations (see, for example, [IH1 [111112011211122]) P rov ide some evidence 
of the existence of a finite-temperature transition for small values of D. On the other 
hand, for large D, even the presence of a finite-temperature transition is not supported 
by numerical simulations [22J. Experiments support the absence of a ferromagnetic 
phase for amorphous systems in general, and provide evidence of a generic glassy 
behavior at sufficiently low temperature. The nature of the transition and of the low- 
temperature phase is however unclear. In particular, no evidence of QLRO has been 
reported. 

The above-reported arguments apply to the RAM in which anisotropy has a 
generic cubic-symmetric distribution. However, in some particular cases, when disorder 
preserves the reflection symmetry s x>a — > —s x>a , s x ^ — > s x> t, for b ^ a (this is, for 
example, realized when the probability distribution vanishes outside the lattice axes), 
there is a standard order-disorder transition with a low-temperature magnetized phase. 
Moreover, a RG analysis of the corresponding Landau-Ginzburg- Wilson theory [2)imi2TI] 
shows that continuous transitions belong to the same universality class as that of the 
random-exchange Ising model (REIM).§ 

In this paper we investigate the critical behavior of the RAM for a uniform 
distribution in the limit D/J — > oo. In this case we can write s x = a x u x , with a x = ±1. 
Thus, the RAM reduces to a particular Ising spin-glass model with Hamiltonian [15] 

it / j Jxy&x&yi Jxy ^x ' ^yi V'J/ 

(xy) 

which we call strong random-anisotropy model (SRAM) (We set J = 1 without loss 
of generality). Model (J5J) differs from the usual ISGM in the bond distribution. Here 

lfil<i>D j g suppressed with respect to the first two analytic terms D 2 and D 4 . 

§ The main difference with respect to the REIM critical behavior (see, e.g., 01|^]) is the approach to 
the asymptotic critical behavior, which is controlled by scaling corrections with exponent A = — a r , 
where a r ~ —0.05 is the specific-heat exponent of the REIM 0. This is much smaller than the 
correction-to-scaling exponent of the REIM, which is Areim ~ 0.25. 
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the random variables j xy on different lattice links are correlated. For instance, one has 
Iln^s/ = 1/27, where the product is over the links belonging to a given plaquette and 
the average is taken with respect to the distribution of the vectors u x . Thus, the SRAM 
is not only of interest as a model of a class of magnetic systems, but also per se, to 
understand how glassy behavior in Ising systems depends on the disorder distribution. 

We study the critical behavior of the SRAM by means of Monte Carlo (MC) 
simulations. Since the model is essentially a spin glass we shall focus on the so-called 
overlap variables ^H] o~ x t x , where a x and t x are associated with two different replicas 
of the model with the same bond variables. For the SRAM one can also consider the 
standard magnetic variables s x = a x u x . The MC results of [22] suggest that these 
quantities are not critical. We confirm here these conclusions so that little will be 
said about magnetic variables in the present study. We will study the behavior of the 
SRAM in the high-temperature phase. This reduces the algorithmic problems — the 
MC algorithm becomes very slow as temperature is reduced — and allows us to consider 
lattices of size L 3 up to L = 30. In order to take into account finite-size effects we use 
the iterative method based on finite-size scaling (FSS) introduced in [22] and applied 
to disordered systems in f2E \ I27 [ 12% } I29j. It allows us to obtain infinite- volume results 
up to £oo ~ 20 (£oo is the infinite-volume second-moment correlation length associated 
with the overlap two-point correlation function) in the high-temperature phase. 

Our MC results do not show a direct evidence of a phase transition for /3 ^ 1.00, 
which is the range of values of /3 that we can reliably simulate. On the other hand, 
our data clearly show FSS as (3 increases, providing indirect evidence of the presence 
of a critical point. Fits of the infinite-volume results obtained for /3 ^ 0.95 indicate 
/3 C = 1.08 ± 0.04. The corresponding critical exponents are: 

Vo = -0.24(4) v = 2.4(6). (4) 

They are defined by x ~ t; 2 ~ Vo and £ ~ (/3 C — /3)~ Uo ; the suffix o is introduced to remind 
that all exponents refer to the overlap variables and not to the magnetic ones. In the 
analysis it is crucial to include corrections to scaling in the FSS (corrections behave as 
L~ w ) and in fits of infinite-volume quantities (nonanalytic corrections behave as £^ as 
£oo — > oo). The exponent u has been determined by studying the critical behavior of 
two universal ratios that involve the four-point and the two-point correlation function 
of the overlap variables. We obtain 

w = 1.0 ±0.4. (5) 

Estimates (J3J) are reasonably close to those obtained for the ISGM (see Table 1 in [30 j 
for a list of recent results) and thus support the conjecture that the SRAM transition 
is in the same universality class as that of the ISGM. Some additional arguments will 
be presented in Sec. El Similar conclusions were reached in [HI] for the two-dimensional 
case. By using the large-cell RG method, it was shown that the two-dimensional SRAM 
has a zero-temperature transition with critical exponents compatible with those of the 
two-dimensional ISGM. 
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The paper is organized as follows. In Sec. 121 we define the quantities that will be 
studied numerically. In Sec.|3]we discuss the MC algorithm. In Sec.EJwe discuss our MC 
simulations, providing evidence for the existence of a critical transition and computing 
the corresponding critical exponents. Finally, in Sec. 03 we draw our conclusions. 

2. Definitions 

We define here the quantities that have been determined in the MC simulation. The 
energy density and specific heat are defined as 

e = Ipwy, C = i/3 2 (W) - W) , (6) 

where V = L 3 is the volume. In our numerical work we focus on the critical behavior 
of the overlap parameter 

q x = a x T Xl (7) 

where a x and r x are two independent replicas of the system with the same couplings 
jxy We consider the correlation function G(x) = (qoq x ), its Fourier transform G(p), the 
corresponding susceptibility x> an d the second-moment correlation length £: 

V 4sin 2 (p min /2) G(p) V } 

where p = (p min , 0, 0), and p min = 2n/L. 

Moreover, we consider quartic correlations of the overlap parameter q x at zero 
momentum. Setting fik = ((^2 x Qx) k ), we define the quartic susceptibilities 

Vxi = ]M - 3^|, Vx22 = tA-Wi 2 ^ (9) 

and the quartic cumulant 

B q = §- 2 . (10) 

Finally, we define the zero-momentum four-point couplings 

Gi " ^ ° 22 " "^' (11) 

We shall also briefly consider magnetic variables associated with s x = o~ x u x . In 
particular, if Uk = (£ si ■ s y ) k ^ 2 ), we define the usual Binder cumulant 

B m = ^ . (12) 

V2 

3. The algorithm 

We consider model Q on simple cubic lattices L 3 with periodic boundary conditions. 
No efficient algorithm is known for generic spin-glass systems (but recently progress 
has been made in some specific cases, see JS3E3E1J and references therein), and thus 
we simply used the Metropolis algorithm with lexicographic choice of the lattice site. 
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In order to reduce the thermalization times, which are the main source of bias in our 
simulations (and also the autocorrelation times, though this is not crucial, since the error 
on the results is mainly due to sample-to-sample fluctuations) we have combined it with 
the random-exchange method — also called parallel tempering or multiple Markov chain 
method — introduced in |3Sl EH EH! (for a recent review with many different applications 
see jHHl; some improvements of the method are discussed in [3§]). 

In practice, we work as follows. We consider Np configurations at different inverse 
temperatures fa (/? = 1/T), i = l,...Np, in a given range [/3 m ; n , /3 max ], that evolve 
according to Hamiltonian (J3J) with the same couplings j xy . Every N ex iterations we 
perform Np — 1 random-exchange moves, trying to swap sequentially fa with fa, then fa 
with fa, up to Pnb-i an d Pnb- Each time we propose a swap of adjacent temperatures 
fa and fa+i, which is accepted with probability Min{l,exp[(^ — E i+ i)(fa — fa+i)]}, 
where E% is the energy of the configuration initially at inverse temperature fa. For each 
sample we run T run Metropolis sweeps on each configuration and then repeat the same 
procedure for iV sam pi e different bond values. Note that, since we need to compute the 
overlap parameter, we simulated at the same time two different replicas of the system. 

Thermalization represents the main source of bias in simulations of disordered 
systems. In all cases we start from a random infinite-temperature spin configuration. 
To check for equilibration, we used the following method. Let Xi be the estimate of the 
overlap susceptibility Xi & t iteration i at the largest (3 value, /3 max , averaged over the 
-Sample disorder realizations. Then, define a block length T Wock and the block-averaged 
quantities 

-I Tbiock 



Xb,; 



Tbiock . 



^2 Xj+(i-l)T hlock ( 13 ) 



In our runs we typically choose Tbi oc k = 5000 or 10000 (all results presented below 
are always in units of Metropolis sweeps). In the tests we report below, in order to 
determine the thermalization times more precisely, we use Tbi oc k = 2000. Finally, plot 
Xb,i as a function of i. For i large, Xb,i becomes constant within error bars, signalling 
equilibration. Data outside the approximately flat region are discarded in the calculation 
of the mean values. 

In the algorithm there are several parameters that must be tuned: /3 m j n , /3 ma x, 
Np, and N cx . Moreover, T run should be large enough to reach equilibration. In our 
simulations the difference A/3 between adjacent (3 values is kept constant, so that 
Anax - Anin = (A^3 — 1) A/?. The highest-temperature value /3 min must be chosen such that, 
at this value of fa the standard Metropolis algorithm is reasonably efficient. In most of 
our simulations we use /? m j n = 0.81. At this value of (3 (it corresponds to an infinite- 
volume correlation length ^ m 3.7) the Metropolis dynamics is reasonably fast. Then, 
at fixed /? m in and /9 max , we investigated how the thermalization time depends on the two 
parameters Np (or, equivalently, A/3, since we work in a fixed (3 interval) and N ex . The 
parameter A/3 controls the acceptance rate of the random-exchange moves and should 
not be too large, otherwise there are no temperature swaps and the exchange dynamics 
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becomes slow. As we shall see, A/3 is not a parameter that requires detailed tuning, 
since there is a somewhat large interval of Af3 values for which the algorithm is equally 
efficient, i.e. that all correspond to the minimal thermalization time (even though the 
acceptance rates change significantly). In order to perform a quantitative comparison, 
we define the thermalization time T t ^ x % as the time such that \xb,i/Xave — 1| = x/100, 
where i = T thjX %/T h i ock and Xave is the average value of \ a ^ Anax- Since the typical 
errors on \ are ~ 1% we consider x% = 2%. The thermalization time T t ^ x % is of 
course a lower bound on the true thermalization time, but it has the advantage of being 
computationally well-defined and thus it allows quantitative comparisons of the results 
corresponding to different choices of the parameters. Moreover, we define a first round- 
trip time Tfrt (in units of Metropolis sweeps): it corresponds to the time needed by a 
given configuration to go at least once through all temperatures. For each simulation we 
collected 2N l 3N samp \ e first round-trip times. Their distribution is not Gaussian, but it has 
an exponential tail e~ cy . Therefore, instead of the average, we found more informative 
to consider the time Tf rt)90 % such that 90% of the measured 2NpN samp \ e first round- 
trip times is smaller than T irtt90 %. In our tests we found Tf rti90 % to be related to the 
thermalization time. Equilibration is reached after a few Tfrt,9o%- 

In Fig.^ for several values of A/?, we show some data corresponding to simulations 
for L = 16, iV sample = 1000, f3 min = 0.81, /3 max = 1.09, and N ex = 10. At the top 
we show the acceptance range. Since A/3 is kept fixed, the acceptance rate is not 
constant. The minimum (maximum) value corresponds to swaps of the configurations 
with the smallest (largest) /3. The acceptance rate for the intermediate (3 values varies 
approximately linearly. In the figure we also show the first round-trip time Tf Ttj90 %, 
the thermalization time T t h,2% (in all cases we used Xave = 644(3) as estimate of % a ^ 
Anax = 1-09), T comp « NpT th)2 %, which is roughly proportional to the computer time. 
The uncertainty on these quantities should be approximately 10%. The data show a 
significant correlation between the thermalization and round-trip times. Their ratio is 
approximately 1.5-2, indicating that a few round trips are needed (and sufficient) to 
thermalize the system. The thermalization time has a minimum in a relatively large 
region of A/3 values, i.e., for 0.02 < A/? < 0.05, corresponding to average acceptance 
rates 10-50%. This suggests that the thermalization time T t h does not depend much 
on the acceptance rate as long as it is not too small (say, larger than 10%). If we 
consider the computer time — and this is what we are really interested in — the optimal 
region is shifted to larger values of A/3. The minimum corresponds to approximately 
A/9 = 0.04, where the acceptance rate varies in the interval [0.115,0.312]. However, an 
additional increase in A/3 does not worsen much the efficiency of the algorithm that is 
nearly optimal even for A/3 = 0.07, where the acceptance rate is rather small (it varies 
between 0.7% and 7%). 

The second interesting parameter is N ex . We compare here two simulations with 
L = 16, iV sam p le = 1000 using N ex = 10, 20 and keeping fixed the other parameters at 
/3 min = 0.81, /3 max = 1.09, and Af3 = 0.04. We find T th , 2 % « 16000 and T frti9 o% « 9000 
for A" cx = 10, and T th2% « 28000 and T frtj90% « 18000 for A" ex = 20. Clearly, N ex = 10 
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Figure 1. We show the acceptance range (above), the quantity iV^Ith 2% (i n units 
of 1000 sweeps) which is roughly proportional to the computer time (middle), the 
thermalization time T tll2 % an d the first round-trip time Tf rt g %, both in units of 1000 
sweeps, (below) versus A/3 for random-exchange runs with L = 16, /3 m i n = 0.81, 
Anax = 1-09, N cx — 10. The numbers reported close to the estimates of T th2 y 
correspond to the number Np of j3 values. 



is better than N ex = 20. This indicates that N ex should be taken relatively small. In 
our simulations we fixed N ex = 10. 

As expected, the thermalization time depends very strongly on /3 max . For example, 
for L = 16, Anin = 0.81, and Af3 = 0.01, we find T thi2% w 16000 for /? max = 1.00 
and T thj2 % ~ 28000 for /3 max = 1.09. Note that that the acceptance rates are similar 
varying from 0.686 (for /3 = 0.81,0.82), to 0.772 (for (5 = 0.99,1.00) and 0.805 (for 
(3 = 1.08,1.09). We have performed an analogous test for L = 24. For j3 min = 0.81, 
A/3 = 0.01, N cx = 10, we find T frti90% w 18000 and T thj2% < 70000 for /3 max = 1.0, 
Tf rt , 90 % w 40000 and T th>2% > 100000 for /? max = 1.05. 

Finally, it is interesting to compare the thermalization times for a Metropolis 
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Table 1. Our runs using the random-exchange method. Here T lun is the number of 
Metropolis sweeps per sample and Tdi sc is the number of discarded Metropolis sweeps. 
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Figure 2. The quartic cumulant B q of the overlap parameter for several lattice sizes. 



simulation without random-exchange moves and a random-exchange simulation. For 
L — 16 and (3 = 1.0 we performed a long Metropolis simulation, considering 5000 
disorder realizations. The thermalization time T t ^ 2 % is approximately 400 • 10 3 . This 
should be compared with a random-exchange run with (3 G [0.81, 1.09], A/3 = 0.04. In 
this case T tn ,2% ~ 16 • 10 3 . Even if the random-exchange run has a larger /3 max , the 
thermalization time is much smaller. 
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Figure 4. The specific heat for several lattice sizes. 



4. Results 

We consider model (J3J) on simple cubic lattices L 3 with periodic boundary conditions. 
We have performed a series of simulations for lattices with 8 < L < 30, using the 
random-exchange method, as explained in Sec. E3 For each value of L we considered a 
range of (3 values such that the thermalization time was less than 10 5 Metropolis sweeps 
per sample. We could therefore take T run « 10 5 . We required this condition in order to 
be able to have N SBmp i e > 1000, and thus precise estimates. Of course, this limits the 
parameter /3 max and, for L = 30, we were able to collect data only up to (3 = 0.95. The 
parameters of our random-exchange runs are reported in Table [0 We are not reporting 
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several other runs with L = 16, iV sam pi e = 1000, that have been used in the numerical 
tests reported in Sec. El Few runs with L = 8, 12, 16 were performed by using only 
Metropolis updatings without random-exchange moves. In Table [Owe report the length 
T run of the run for each value of the parameters and T^[ SC , the number of sweeps that 
have been discarded before measuring (this parameter has been chosen conservatively 
to avoid any thermalization bias). Simulations took approximately 2.5 CPU years of a 
workstation equipped with an AMD Opteron Processor 246 (2 GHz clock). 

In Figs. El and El we show the quartic cumulant B q defined in Eq. (TTTTf) and the 
ratio £/L for several lattice sizes L. In the region of L and (3 covered by our MC data, 
neither B q nor £/L have a crossing point. Thus, we do not have direct evidence for a 
finite-temperature transition in the range of temperatures that we can reliably simulate. 
No indication of a phase transition in this range of /3 is also provided by the results for 
the specific heat — they are shown in Fig. 0] — and by magnetic variables — there is no 
indication of a crossing point in the results for B m defined in Eq. (fT2"|) . see Fig. 

These results are consistent with two possible scenarios: (i) the system becomes 
critical at j3 = (3 C with f3 c ^ 1.05 (with the possibility of a zero-temperature transition, 
i.e. p c = +oo); (ii) the system never shows criticality and even for (3 = oo the 
correlation length is finite. We will now show that our data allow us to exclude this 
second possibility since, as j3 increases, the model shows critical behavior; more precisely, 
our data for the overlap variables show FSS as expected close to a critical point. In 
order to make the check as reliable as possible we will study the FSS behavior of the 
ratios A({3, sL)/A((3, L), where A(/3, L) is a long-distance quantity and s is an arbitrary 
number. In the FSS limit (i.e. for L,^(j3,L) — > oo at £(/3,L)/L fixed) we should have 

A{f3,sL) 



A(/3,L) 



F A [ 8 ,Z(J3,L)/L] 



(14) 
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Figure 6. The FSS curve of the susceptibility x for s = 3/2. 



where Fa(s, z) is a universal function. Note that in this formulation there are no free 
parameters and thus one can make an unbiased test of FSS. 

The FSS curves computed by using Eq. (|T3|) can be used to determine infinite- 
volume quantities for large values of the correlation length. For this purpose we employ 
the extrapolation method of J2E] (see also JH3 E] for a discussion of the efficiency of 
this technique). Indeed, in the absence of scaling corrections, Eq. (|H|) allows us to 
compute A(/3, sL) on a lattice of size sL in terms of quantities defined on a lattice 
of size L and of the function Fa(s,z). In practice, one works as follows. First, one 
performs several runs, determining A(/3,sL), A(/3,L), £(/3,sL), and £(/3, L). By means 
of a suitable interpolation, this provides the FSS function Fa(s, z) for A and £. Then, 
Aoo(f3) and £oo(/3) are obtained from A(P,L) and £(/?, L) by iterating Eq. (|H|) and the 
corresponding equation for £(/3, L). Of course, one must be very careful about scaling 
corrections, discarding systematically lattices with small values of L till results become 
independent of L within error bars. 

In Figs. El El and El we plot the ratios ^U^i for \i £> G<±, and G22, fixing s = 3/2. 
The curves for x and £ show significant scaling corrections. However, they apparently 
decrease very rapidly with L and indeed the data corresponding to the pairs L = 16, 24 
are only slightly different from those with L = 20, 30. Therefore, they strongly suggest 
that FSS holds, although, for L ^ 30, scaling corrections are significant compared to 
our error bars (for x our data have a relative error of approximately 1% for L = 20 
and L = 30). For G4 and G22 corrections are apparently weaker and indeed the 
data corresponding to L — 20, 30 are compatible with those with L = 12, 18 and with 
L = 16, 24. In this case FSS holds within error bars. 

The presence of scaling corrections in the FSS curve for £ does not allow us to use 
straightforwardly the iteration method of [25/ and makes it necessary to include scaling 
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Figure 7. The FSS curve of the correlation length £ for s = 3/2. 
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Figure 8. The FSS curve of the zero-momentum quartic coupling G4 for s = 3/2. 



corrections in the scaling Ansatz. As done in [27] , we use a more general Ansatz of the 
form 

A(J3,sL) 



F A (s } £(/?, L)/L) + L-"G A (s, £(J3, L)/L), 



(15) 



A(P,L) 

where Ga{s, z) is a new universal scaling function and u is a to-be-determined correction- 
to-scaling exponent. 

Our data are not precise enough to allow a determination of uj. However, it is easy 
to convince oneself that uo cannot be arbitrarily large. If we consider the correlation 
length £(/3, L), the correction-to-scaling exponent should be less than 2, since corrections 
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Figure 9. The FSS curve of the zero-momentum quartic coupling G22 for s — 3/2. 
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with exponent 2 appear necessarily j42j. They are related to the very definition of the 
second-moment correlation length we use on a finite lattice. Indeed, it corresponds to 
its infinite-volume counterpart only up to terms of order L~ 2 . For x, G4, and G22 there 
are corrections due to the analytic background |l2l IS] that are proportional to L 2 ~ Vo , 
where r\ is the susceptibility exponent, x ~ £ 2 ~ Vo - As we shall see, r\ w —0.3, so 
that c<j < 2.3. On the other hand, one cannot set a priori a lower bound on uj and, in 
principle, uj can be arbitrarily small. We will present below analyses with u as small 
as 0.2, presenting results for to = 0.2, u — 1, and to — 2. This choice will be justified 
below. 
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Figure 11. Log- log plot of the infinite- volume estimates of %oo versus £ 



In order to determine r\ , we consider the data reported in Figs. El and [7| and 
apply the extrapolation technique using Eq. (fT5|) . As in J2HJ the scaling curves are 
parametrized as polynomials in exp[— £(/3, L)/L] satisfying Ga(s,0) = ^4(5,0) = 1. A 
fourth-order polynomial for F^(s,z) and a third-order polynomial for Ga{s,z) appear 
to be adequate. In Fig. HO we show the infinite- volume correlation length ^ (/?). It can 
be determined reliably up to (3 = 0.95, essentially because we have results for the largest 
lattice, L = 30, only up to this value of [3. The dependence on uj is not large, but in 
any case larger that the statistical errors. For j3 = 0.95, we find £<*, = 26.4(1.7), 21.5(5), 
19.4(4) for uj = 0.2, 1, 2, respectively. Whatever is the value of uj, £oo(/3) increases quite 
rapidly, confirming that the system eventually becomes critical. Evidence of criticality 
is also provided in Fig. ITT1 where we show a log-log plot of Xoo versus £oo. The behavior 
is linear, indicating that 

xoo ~ &-* (i6) 

where r\ is a critical exponent. Note that the dependence on uj of Xoo at fixed ^ is 
much smaller than the uj dependence of Xoo an d £00 at fixed f3. This is due to the fact 
that, at fixed /?, Xoo and £oo are strongly correlated: they both increase with decreasing 
uj, in such a way that the ratio Xoo/£,^ Vo has a tiny dependence on uj. To obtain an 
estimate of r\ we fitted Xoo and ^ to lnx^ = a + (2 — ^ln^x,. To estimate the 
scaling corrections, the fit has been repeated several times, each time including only the 
data satisfying (3 > /? m i n - Results are reported in Table 121 for different values of uj. In 
order to have an additional check on the stability of the results, we have also repeated 
the analysis twice: we present results obtained by using all data (L min = 8) and results 
obtained using only lattices with L > 12 (L min = 12). 

The results presented in Table El depend somewhat on uj and /3 m i n . At fixed uj, r\ 
decreases with increasing /3 m j n , reaching an approximate plateau within error bars for 
Anm ~ 0.84. At fixed /? m i n the estimates decrease with increasing uj. A conservative 
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Table 2. Estimates of the exponent r\ Q for several values of w and L Tl 







-^min O 






-^min *■& 




Pmin 


a; = 0.2 


LU= 1 


uo = 2 


cj = 0.2 


U= 1 


u = 2 


0.76 


-0.161(7) 


-0.177(6) 


-0.185(6) 


-0.148(9) 


-0.159(8) 


-0.167(7) 


0.78 


-0.191(9) 


-0.203(8) 


-0.214(7) 


-0.175(11) 


-0.185(9) 


-0.196(9) 


0.80 


-0.210(12) 


-0.221(10) 


-0.233(9) 


-0.196(15) 


-0.203(12) 


-0.216(11) 


0.82 


-0.221(15) 


-0.234(12) 


-0.247(10) 


-0.205(21) 


-0.214(15) 


-0.227(14) 


0.84 


-0.236(22) 


-0.251(16) 


-0.268(14) 


-0.225(31) 


-0.234(22) 


-0.242(19) 


0.86 


-0.233(32) 


-0.258(22) 


-0.276(19) 


-0.217(50) 


-0.233(31) 


-0.251(27) 


0.88 


-0.236(50) 


-0.265(31) 


-0.288(26) 


-0.213(81) 


-0.241(45) 


-0.262(38) 
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Figure 12. Estimates of Ga,oo versus £oo for three different values of u>. 



estimate can be obtained by noting that all results with 0.84 < p min < 0.86 — for larger 
values of (3 m \ n error bars are quite large — lie in the interval —0.28 ^ r\ < —0.18, 
including statistical errors. We thus end up with the following estimate: 

Vo = -0.23(5). (17) 

Then, we consider the RG-invariant quantities G4 and G^- Again, we use Eq. (J15|) and 
the iteration algorithm to determine infinite- volume estimates G± t00 {f3) and (^22,00 (/^)- 
Reasonable results are only obtained if the data with L = 8 are discarded, i.e. if we 
only consider L > L m [ n = 12. Indeed, if all MC data are used in the FSS extrapolation, 
G^oo(P) and G22,oo (P) apparently do not converge to a finite value as (3 — > cxd. For 
-^min = 12, the u dependence is of the order of the statistical error bars: for instance, 
for (3 = 0.95 we predict G 4 ,oo = 85(10), 91(7), 95(4), and G 22 ,oo = -13.4(7), -13.1(7), 
— 11(2), for u = 0.2, 1, 2. The results are shown in Figs. IT2l and fT3l They show a 
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rapid variation with ^ for small values of the correlation length and reach a plateau 
for £00 > 15. In order to estimate the critical value for ^ — ■> 00, we perform a fit of the 
form 



ft„ 



ft* + 



■n ' 

00 



;is) 



where 7?.*, a^, and Q are free parameters, TZ = G4 or G22- Results are reported in Tables 
El and |U The dependence on u is here quite small, as it has to be expected since there 
is no evidence of scaling corrections in the FSS curves for G\ and Gw The dependence 
we observe here is mainly related to the corrections affecting the correlation length that 
is also used in the extrapolation. Again our final estimate is obtained conservatively by 
looking at the variation with u and f3 n 



For Gl, we take as our final estimate 



G* 4 = 88(8), 



(19) 



where the error is such to include all results with /? min > 0.81. Analogously, we estimate 

G* 22 = -11(4). (20) 

The fit parameter Q that appears in Eq. (JT%|l is not completely arbitrary and can be 
related to the exponent u used in the FSS extrapolation. Thus, the results of Tables El 
and |U can be used to obtain constraints on the value of u. In general, given an infinite- 
volume RG-invariant quantity 7Z{/3), close to a critical point we expect an expansion of 
the form: 

TZ(P)=n*(l + a\P-p c \ A ), (21) 

with a positive correction-to-scaling exponent A. According to the RG theory, scaling 
corrections may have several origins |43| l4l)j : 

(i) There are analytic corrections of the form \/3 — /3 c | n , n being an integer. 
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Table 3. Estimates of G\ and for different values of ui. Here L m ; n = 12. 





u! 


= 0.2 


U) 


= 1 


U) 


= 2 


Pmin 


(_T 4 


n 


(_T 4 


Q 


(jr 4 


Q 


0.780 


58(26) 


0.36(19) 


72(13) 


0.47(17) 


78(8) 


0.55(15) 


0.790 


74(13) 


0.58(23) 


78(9) 


0.58(19) 


83(6) 


0.70(18) 


0.800 


85(7) 


0.95(30) 


86(6) 


0.89(25) 


87(4) 


0.93(21) 


0.810 


87(6) 


1.08(34) 


87(5) 


1.01(29) 


88(4) 


0.98(23) 


0.815 


86(7) 


1.00(38) 


88(5) 


1.07(34) 


88(4) 


0.98(27) 


0.820 


89(6) 


1.25(45) 


89(5) 


1.20(38) 


89(4) 


1.07(29) 


0.825 


88(7) 


1.13(50) 


88(6) 


1.06(40) 


88(5) 


1.01(32) 


0.830 


91(6) 


1.47(60) 


90(5) 


1.34(49) 


90(4) 


1.20(38) 


0.835 


91(6) 


1.49(76) 


89(6) 


1.18(56) 


89(5) 


1.07(43) 


0.840 


95(4) 


2.49(1.25) 


92(5) 


1.66(78) 


91(4) 


1.37(54) 


0.845 


89(10) 


1.13(95) 


86(11) 


0.82(65) 


84(11) 


0.70(51) 


0.850 


91(8) 


1.45(1.20) 


87(9) 


0.94(73) 


85(10) 


0.76(55) 



111 



There are nonanalytic corrections related to the irrelevant operators. They have 
the form \(3 - p c \ nuuJ \ \/3 - p e \™->i*™*>i , etc., where u t = -y t and y t are the RG 
dimensions of the irrelevant RG operators, and n, m are integers. 

There are corrections related to the analytic background. For instance, in the 
susceptibility there are corrections proportional to |/3 — /3 C | 7 . 

In jlTj an argument was given to show that renormalized coupling constants, like G 4 
and G22, do not have analytic corrections. The absence of this type of corrections was 
verified in the continuum O(N) 4 theory in d dimensions to order 1/N jUj and in the 
two-dimensional Ising model ^HllMl- I n the absence of analytic corrections, the leading 
contribution is due to the irrelevant operators, and therefore A = oj\v, where uo\ = —y\ 
and ?/i is the RG dimension of the leading irrelevant operator in the model. Therefore, 
in Eq. ()18j) the exponent Q should be identified with uj\. 

Corrections appearing in Eq. (|2T|) are strictly related to corrections appearing in 
FSS. In the FSS case no analytic corrections are expected, not only in G 4 or G 2 2 but in 
any quantity. || Corrections of type (ii) correspond to terms of the form L~ nuJi , ^- nw »- maJ j ( 
so that the leading correction in FSS has the form L~ Wl . Thus, the exponent Q should 
be identified with the exponent that controls FSS corrections and has been used in the 
FSS Ansatz (|T3jl. The results reported in Tables El and |U give estimates of Q that show 
a tiny dependence on to and we can safely estimate 



Q = 1.0(4). 



(22) 



|| A careful discussion is presented in |21|. The absence of analytic corrections in FSS is shown in j^D] 
and has been checked to order L~ 2 in the two-dimensional Ising model jj^l an d to order L~ l in the 
three-dimensional XY model 51 . 
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Table 4. Estimates of G\^ and fl for different values of u. Here L m - m = 12. 



uJ 



0.2 



u> 



UJ 



22 



n 



22 



n 



22 



n 



0.780 


-8.1(3.8) 


0.69(22) 


-10.6(1.8; 


0.81(16) 


-10.1(1.5; 


0.76(13) 


0.790 


-10.2(3.0) 


0.88(27) 


-11.5(1.6; 


0.94(19) 


-10.9(1.3; 


0.88(16) 


0.800 


-11.1(2.8) 


1.00(32) 


-11.8(1.5; 


1.01(22) 


-11.3(1.3; 


0.94(18) 


0.810 


-10.6(3.3) 


0.92(34) 


-11.7(1.6; 


0.99(24) 


-11.1(1.4; 


0.90(19) 


0.815 


-10.1(3.8) 


0.86(37) 


-11.8(1.7; 


1.01(27) 


-10.9(1.6; 


0.87(21) 


0.820 


-11.6(2.9) 


1.10(44) 


-12.3(1.5; 


1.14(29) 


-11.6(1.3; 


1.01(24) 


0.825 


-9.1(5.4) 


0.73(44) 


-11.5(2.0; 


0.93(30) 


-10.2(2.1; 


0.75(24) 


0.830 


-11.2(3.6) 


1.03(54) 


-12.4(1.6; 


1.18(37) 


-10.9(1.8 


0.87(27) 


0.835 


-9.6(5.6) 


0.78(55) 


-11.5(2.3; 


0.91(38) 


-10.0(2.7 


0.71(30) 


0.840 


-12.8(2.9) 


1.53(93) 


-12.6(1.7; 


1.25(49) 


-11.1(2.0; 


0.90(35) 


0.845 


-9.5(7.2) 


0.73(71) 


-10.2(4.0; 


0.67(45) 


-7.9(5.7) 


0.49(36) 


0.850 


-10.0(6.6) 


0.81(80) 


-10.3(4.0; 


0.68(47) 


-8.1(5.7) 


0.51(38) 



This result justifies a posteriori our choice of not considering values of uj smaller than 
0.2. The estimate of Q restricts the interval in which uj can vary. This allows us to 
obtain more precise estimates for n and G% 2 (the estimate of G\ does not vary since 
this quantity has a very small dependence on to). By assuming 0.6 < uj < 1.4 we obtain 
the results: 

Vo = -0.24(4), (23) 

G* 22 = -11.5(2.5). (24) 

Finally, we try to determine the position of the critical point. For this purpose, we 
analyze Xoo and £oo separately as 

logXoc(/3) =A X - 7olog(/3 c - P), 

log £00 03) =As-v \og({3 c -{3), (25) 

where A x , A%, j , u , and (3 C are free parameters. As before, we repeat the fit several 
times, including each time only the data with (3 > (3 m \n ■ Results using the infinite- 
volume data obtained with u = 1 are reported in Table El The estimates of j3 c show a 
tiny dependence on /3 m i n ; moreover, the results obtained by using Xoo and £00 are fully 
consistent. We take as our final estimate 

P c = 1.08 ±0.01 ±0.03, (26) 

where the first error is statistical and the second one is systematic and takes into account 
the change in the estimate as u varies between 0.6 and 1.4. The data indicate therefore a 
transition point that is just beyond the interval of /3 where we have made the simulations 
and thus they are compatible with the absence of a crossing point in B q and £/L and with 
the fact that the estimates of these two quantities at fixed L get closer as (3 increases. 
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Table 5. Estimates of /3 C and 7 (second and third column) obtained from the fit of x 
and of P c and v (fourth and fifth column) obtained from the fit of £oo . The infinite- 
volume estimates of Xoo and £oo have been obtained by applying the extrapolation 
method based on Eq. I|15|l with u> = 1 to all data (L m [ n = 8). 





fit of 


Xoo 


fit of £oo 


Pmin 


& 


lo 


ft 


Vo 


0.76 


1.060(2) 


4.78(4) 


1.083(3) 


2.43(4) 


0.78 


1.059(2) 


4.76(5) 


1.075(4) 


2.33(4) 


0.80 


1.060(2) 


4.76(6) 


1.073(4) 


2.29(5) 


0.82 


1.062(3) 


4.84(8) 


1.073(5) 


2.30(6) 


0.84 


1.064(4) 


4.89(12) 


1.072(6) 


2.28(8) 


0.86 


1.070(6) 


5.11(19) 


1.076(9) 


2.35(13) 


0.88 


1.074(9) 


5.27(32) 


1.080(13) 


2.42(21) 


0.90 


1.082(16) 


5.60(62) 


1.080(19) 


2.41(34) 



We can also perform a more quantitative check by using the results of [30J for the 
critical-point values B* and (£/L)*. They quote: B* = 1.475(6) (bimodal distribution) 
and B* = 1.480(14) (Gaussian distribution); (£/£)* = 0.627(4) (bimodal distribution) 
and (£/£)* = 0.635(9) (Gaussian distribution). These results are compatible with ours 
for B q and £/L close to (3 = 1.08. For (3 = 1.07 we have B q = 1.411(4) (L = 12), 
B q = 1.434(6) (L = 16), and £/L = 0.662(4) (L = 12), £/L = 0.648(4) (L = 16). 
These results are very close to the ISGM estimates and show the correct trend. They 
are therefore consistent with the conjecture that the ISGM and the SRAM belong to 
the same universality class. Magnetic variables behave differently: as is clear in Fig. 
the estimates at fixed L are well separated, so that we do not expect any crossing at 
j3 ~ 1.08. Thus, no criticality is expected in magnetic quantities: below the transition 
temperature the system is still paramagnetic. 

Using the results reported in Table we can also estimate the critical exponents 
7o and v . We obtain: 

u = 2.4 ± 0.2 ± 0.4 (27) 

7 = 5.3 ±0.3 ±1.0. (28) 

The estimates correspond to /5 min = 0.88. The first reported error is the statistical error 
and should be large enough to include most of the dependence on /3 m i n . The second one 
gives the variation of the estimate with u. 

5. Conclusions 

In this paper we study the SRAM by MC simulations with the purpose of clarifying 
whether this model shows a critical behavior. We find clear evidence of FSS and 
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criticality, with a critical point located at 

P c = 1.08(4). (29) 

Correspondingly, we compute the critical exponents associated with the critical behavior 
of overlap observables. We obtain: 

Vo = -0.24(4) (30) 

v = 2.4(6) (31) 

7o = 5.3(1.3). (32) 

These estimates are close to those obtained for the ISGM. Ref. J2H] quotes v = 2.22(15), 
r] = -0.349(18), while [30! reports u = 2.39(5) and rj = -0.395(17). Other estimates 
are reported in Table 1 of [201 • With the quoted error bars there is a small discrepancy 
between our estimate of r\ and those of [2BJ ED] • This difference should not be taken too 
seriously, since there are similar discrepancies among the estimates obtained by different 
groups for the bimodal ISGM, see Table 1 of [HOj. For instance, our result is close to 
rj = —0.26(4) reported in [26, • Note that J2B| also reports an estimate of the subleading 
exponent u = 0.7(3). This result is obtained from fits of Xoo vs £oo- in this case, 
however, analytic corrections are expected and therefore the leading scaling correction 
is £i> ~ Co°' 45 - Therefore, the result obtained in J2H] provides only a strong indication 
that to > \jv but does not give a quantitative estimate of the exponent associated with 
the leading irrelevant operator. The bound uj > \/u is fully consistent with our results. 

The possibility that the SRAM is in the same universality class of the ISGM was 
put forward in [14] and found to be consistent with numerical results in two dimensions 
in J3T]. It looks very plausible, since the SRAM is nothing but an Ising model with local 
disorder and frustration. In some sense we can think of the SRAM and of the ISGM as 
two different versions of the same model: in the SRAM disorder is associated with lattice 
sites, while in the ISGM disorder is associated with lattice bonds. They are analogous 
to the site-diluted and bond-diluted Ising model, whose Hamiltonian is given by Eq. (JHJ) 
with j xy = r x r y (site dilution) and j xy = r xy (bond dilution), r being a random variable 
such that r = 1 (r = 0) with probability p (resp. 1 — p). Note that the site-diluted 
model can also be seen as a bond-diluted model with a correlated bond distribution, 
exactly as is the case for the SRAM. Nonetheless, there is little doubt — though a precise 
numerical check is still missing — that the two models belong to the same universality 
class. 

Note that the SRAM is less frustrated than the standard ISGM, since in the SRAM 
Iln^s/ = 1/27. This fact does not rule out our conjecture since it is known that maximal 
frustration is not necessary to obtain glassy behavior. For instance, the random-bond 
Ising model with j xy = +1 with probability p and j xy = — 1 with probability 1 — p has 
a glassy low-temperature phase for jH2| 0.233 ^ p ^ 0.767. 

It is interesting to generalize to SRAM by considering N- dimensional unit vectors 
Ui. The correlation of the bond variables around a lattice plaquette becomes Iln^y = 
1/iV 3 , which implies that bond correlations vanish for iV — » oo. Thus, for iV = oo, the 
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SRAM is just an ISGM with a different continuous bond distribution. In this limit, 
therefore, the two models belong to the same universality class. If our conjecture for 
N = 3 is correct, the same should hold for any N > 3. For N = 1 it is enough to 
redefine <7j — > UiOi to reobtain the standard ferromagnetic Ising model. The behavior 
for N = 2 is unclear, since this model is less frustrated than that with N = 3 studied 
here. In analogy with what happens in the random-bond Ising model mentioned in the 
previous paragraph, it could have a glassy transition or a ferromagnetic transition (the 
nature of the ferromagnetic transition is still object of debate, see [HH]). 

Our results are not precise enough to show conclusively that the SRAM and the 
ISGM belong to the same universality class. Because of the somewhat large scaling 
corrections, one needs to perform simulations on larger lattices and deeper in the critical 
regime (note that for L = 30, our largest lattice, our data extend only up to (3 = 0.95, 
that is quite far from the critical point) to obtain precise and reliable estimates of 
the critical exponents. In particular, it would be interesting to see whether the small 
difference in the estimates of rj disappears (at the same time it is important to include 
scaling corrections in the analysis of the ISGM results, to understand the reliability of 
the quoted error bars). Since we have reasonably precise estimates for G\ and G\ 2 -, it 
is also interesting to compute the same quantities in the ISGM. This would provide an 
additional check of the conjecture. Work in this direction is in progress. 

A question that remains open is the behavior of the RAM for finite anisotropy 
D. If there is indeed a low-temperature phase with QLRO for small D as predicted 
in fH\ IT2"| IT3*] . then there should be a critical value D* such that ISGM behavior is 
observed only for D > D*. Nothing is known about D* and we cannot even exclude 
that D* = oo, so that ISGM behavior is observed only for model ( |3~j) . 
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